//Michael Alstad
//Extra Lab
//Throwing Rock Exercise
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;

float calcPos(float, float, float, float);
float quadFormula(float, float, float);
float percentError(float, float);
float convertToRad(float);


//currently does not calculate which angle gives the greatest distance
int main()
{
	const float Gravity = -9.80665;
	//input variables
	float fInitHeight;
	float fAngle;
	float fInitSpeed;
	float fTimeIncrement;

	//initial horizontal and vertical speed
	float fInitHSpeed;
	float fInitVSpeed;

	//calculated variables
	float fMaxVert = 0;
	float fVSpeed;
	float fVert;
	float fDist = 0;
	float fTime = 0;

	//theoretical variables
	float fTrueTime;
	float fTrueTimeAtMaxVert;
	float fTrueDist;
	float fTrueMaxVert;

	float fPercentErrorTime;
	float fPercentErrorDist;
	float fPercentErrorMaxVert;

	//get input data
	cout << "Input throwing height in meters: ";
	cin >> fInitHeight;
	cout << "Input throwing angle in degrees: ";
	cin >> fAngle;
	cout << "Input initial speed in m/s: ";
	cin >> fInitSpeed;
	cout << "Input time increment in seconds: ";
	cin >> fTimeIncrement;

	//calculate initial horizontal and vertical speeds
	fInitHSpeed = fInitSpeed  * cos(convertToRad(fAngle));
	fInitVSpeed = fInitSpeed * sin(convertToRad(fAngle));
	//gives variables initial values before entering update loop
	fVSpeed = fInitVSpeed;
	fVert = fInitHeight;

	//fTrueTime represents the time it takes for the rock to fall to the ground
	fTrueTime = quadFormula(Gravity * .5, fInitVSpeed, fInitHeight);
	fTrueTimeAtMaxVert = -1 * fInitVSpeed / Gravity;
	
	//calculates the horizontal distance reached at fTrueTime
	fTrueDist = calcPos(0, fInitHSpeed, fTrueTime, 0);
	fTrueMaxVert = calcPos(fInitHeight, fInitVSpeed, fTrueTimeAtMaxVert, Gravity);

	fTime += fTimeIncrement;

	//begin simulation with current angle		
	for(fTime; fVert > 0; fTime += fTimeIncrement)
	{
		//updates vertical speed
		fVSpeed += Gravity * fTimeIncrement;
		//updates vertical position
		fVert += fVSpeed * fTimeIncrement;
		//updates horizontal position
		fDist += fInitHSpeed * fTimeIncrement;

		if(fVert > fMaxVert)
			fMaxVert = fVert;
		//cout << "Time: " << fTime << " seconds " << endl;
		//cout << "Height: " << fVert << " meters" << endl;
		//cout << "Distance: " << fDist << " meters" << endl << endl;

	}

	//calculate error percentage by comparing theoretical values to the simulated ones
	fPercentErrorTime = percentError(fTrueTime, fTime);
	fPercentErrorDist = percentError(fTrueDist, fDist);
	fPercentErrorMaxVert = percentError(fTrueMaxVert, fMaxVert);
	
	//output results and margin of errors
	cout << setprecision(3) << fixed;
	cout << "\nTheoretical Total Time: " << fTrueTime << endl;
	cout << "Simulated Total Time: " << fTime << endl << endl;

	cout << "Theoretical Distance: " << fTrueDist << endl;
	cout << "Simulated Distance: " << fDist << endl << endl;

	cout << "Theoretical Maximum Height: " << fTrueMaxVert << endl;
	cout << "Simulated Maximum Height: " << fMaxVert << endl << endl;

	cout << setprecision(2) << fixed;
	cout << "\nPercent of Error for time: " << fPercentErrorTime << "%\n";
	cout << "Percent of Error for distance: " << fPercentErrorDist << "%\n";
	cout << "Percent of Error for maximum height: " << fPercentErrorMaxVert << "%\n";

	return 0;
}

//formula to calculate theoretical horizontal and vertical positions
float calcPos(float init, float initVel, float time, float accel)
{
	//f(t) = f(0) + v(0)t + 1/2at^2
	return init + initVel * time + 0.5 * accel * time * time;
}

float quadFormula(float A, float B, float C)
{
	//quadratic formula really has two answers but for this project
	//one answer will be negative and the other should be positive
	float ans1 = (B * -1 + sqrt((B * B) - (4 * A * C))) / (2 * A);
	float ans2 = (B * -1 - sqrt((B * B) - (4 * A * C))) / (2 * A);

	//returns the larger (positive) value
	if(ans1 > ans2)
		return ans1;
	else
		return ans2;
}

//determines percent of error between
float percentError(float theoretical, float simulated)
{
	float answer = 100 * ((theoretical - simulated) / theoretical);

	if(answer < 0)
		answer *= -1;

	return answer;
}

float convertToRad(float angle)
{
	return angle * (3.14159265 / 180);
}
